Population \((N = 300)\)
Transmission/Infection Rate \((r = 0.02)\)
Percentage of healthy individuals with symptoms \((g = 0.2)\)
Initial number of infected individuals day1 \((d_1 = 5)\)
Initial number of cured individuals day1 \((c_1 = 0)\)
Percentage of individuals receiving Treatment A \((PTA = 1)\)
Percentage of individuals receiving Treatment B \((PTB = 0)\)
Treatment A effectiveness \((A_{eff} = 0.75)\)
Treatment B effectiveness \((B_{eff} = 0.25)\)
Treatment A cost \((A_{cost} = 20.00)\)
Treatment B cost \((B_{cost} = 15.00)\)
Number of days since first case \((day = 38)\)
Model upgrade with randomness - Additional Variables
Initial number of known cases (untreated) day1 \((u_1 = 0)\)
Percentage of individual being tested daily \((PS = 1)\)
Test sensitivity \((tp = 0.99)\)
Test specificity \((tn = 0.95)\)
Testing cost \((Test.cost = 1)\)
complex.model.rand = function(N=300, r=0.02, g=0.2, d1=5, c1=0, u1=0, PS=1, tp=0.99, tn=0.95, PTA=1, PTB=0, A.eff=0.75, B.eff=0.25, day = -1, A.cost = 20, B.cost = 15, Test.cost = 1){
# Initial conditions
curr.day = 1
d = d1
c = c1
u = u1
# Calculate initial statistics
h = N - d - c
hs = round(g * h)
s = hs + d
a = s - u
ns = round(PS * a)
ad = round(PS * (d - u))
as = ns - ad
tps = rbinom(1, as, 1-tn)
tpd = rbinom(1, ad, tp)
tpt = tps + tpd
tpt1 = tpt
NDAT = tpd + u
NDTA = 0
NDTB = 0
DCA = 0
DCB = 0
NAT = 0
NTA = 0
NTB = 0
SDC = 0
New = min(N - c - d, round(r * hs * d))
# Table storing S, I, R compartments over time
df = data_frame(Day = curr.day,
Healthy = h,
Symptopms = hs,
Diseased = d,
Tot.Symp = s,
Avail.screen = a,
Untreated = u,
Screen.Tot = ns,
Screen.Sym = as,
Screen.Dis = ad,
Test.Pos.Sym = tps,
Test.Pos.Dis = tpd,
Test.Pos.Total = tpt,
Avail.TRT.Dis = NDAT,
TRT.A.Dis = NDTA,
TRT.B.Dis = NDTB,
TRT.A.Cure = DCA,
TRT.B.Cure = DCB,
Avail.TRT.Tot = NAT,
TRT.A.Tot = NTA,
TRT.B.Tot = NTB,
Sum.Dis.Cured = SDC,
Cured = c,
Catch.Dis = New)
# Predict S, I, R compartments by day
while(day == -1 & d > 0 | (day != -1 & curr.day < day)){
# if number of days is specified, calculate up until that day
# else calculate until there is no more diseased individual
# Update current day and c, u, d
curr.day = curr.day + 1
c = c + SDC
u = d - SDC
d = d + New - SDC
# Calculate other variables
h = N - d - c
hs = round(g * h)
s = hs + d
a = s - u
ns = round(PS * a)
ad = round(PS * (d - u))
as = ns - ad
tps = rbinom(1, as, 1-tn)
tpd = rbinom(1, ad, tp)
tpt = tps + tpd
NDAT = tpd + u
NDTA = round(PTA * NDAT)
NDTB = min(NDAT - NDTA, round(PTB * NDAT))
DCA = min(NDTA, round(NDTA * A.eff))
DCB = min(NDTB, round(NDTB * B.eff))
if (curr.day == 2){
NAT = tpt + tpt1
}
else{
NAT = tpt + u
}
NTA = round(PTA * NAT)
NTB = NAT - NTA
SDC = DCA + DCB
New = min(N - c - d, round(r * hs * d))
data = c(curr.day, h, hs, d, s, a, u,
ns, as, ad, tps, tpd, tpt,
NDAT, NDTA, NDTA, DCA, DCB, NAT, NTA, NTB,
SDC, c, New)
df = rbind(df, data)
}
df$Cost = df$TRT.A.Tot * A.cost + df$TRT.B.Tot * B.cost + df$Screen.Tot * Test.cost
return(df)
}
model.plot = function(model){
cost = sum(model$Cost)
sick = sum(model$Diseased)
peak = max(model$Diseased)
ggplot(data = model, aes(x = Day)) +
geom_line(aes(y = Healthy, color = "Healthy")) +
geom_line(aes(y = Diseased, color = "Diseased")) +
geom_line(aes(y = Cured, color = "Cured")) +
scale_color_discrete(name = "Compartment") +
labs(y = "Population", title = paste("Complex SIR Model: infected = " , sick, ", peak = ", peak, ", cost = $", cost, sep = "")) +
theme_minimal()
}
df = complex.model.rand()
df
model.plot(df)
line.plot = function(model){
list(geom_point(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.8),
geom_point(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.8),
geom_point(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.8),
geom_line(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.1),
geom_line(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.1),
geom_line(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.1))
}
model.rep.plot = function(reps){
plot = ggplot() + labs(y = "Population", title = paste("Complex SIR Model with", reps, "reps")) + theme_minimal()
for(i in 1:reps){
df = complex.model.rand()
plot = plot + line.plot(df)
}
plot = plot + scale_color_discrete(name = "Compartment")
plot
return(plot)
}
model.rep.plot(10)
model.rep.plot(100)
accumulate_by <- function(dat, var) {
var <- lazyeval::f_eval(var, dat)
lvls <- plotly:::getLevels(var)
dats <- lapply(seq_along(lvls), function(x) {
cbind(dat[var %in% lvls[seq(1, x)], ], frame = lvls[[x]])
})
dplyr::bind_rows(dats)
}
visual_data = select(df, c(Day, Healthy, Diseased, Cured)) %>% tidyr::gather(key = Compartment, value = Population, Healthy, Diseased, Cured) %>% accumulate_by(~Day)
plot_ly(data = visual_data, type = "scatter", mode = "lines",
x = ~Day, y = ~Population, color = ~Compartment, frame = ~frame) %>%
layout(hovermode = 'compare') %>%
animation_opts(frame = 100, transition = 0, redraw = FALSE, mode = "immediate") %>%
animation_slider(hide = T) %>%
animation_button(x = 1, xanchor = "right", y = -0.1, yanchor = "bottom", label = "View Plot")
`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
[90mThis warning is displayed once every 8 hours.[39m
[90mCall `lifecycle::last_warnings()` to see where this warning was generated.[39m